Overcoming Gamma Ray Constraints with Annihilating Dark Matter 

in Milky Way Subhalos 
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We reconsider Sommerfeld-enhanced annihilation of dark matter (DM) into leptons to explain 
PAMELA and Fermi electron and positron observations, in light of possible new effects from sub- 
structure. There is strong tension between getting a large enough lepton signal while respecting 
constraints on the fluxes of associated gamma rays; we show how DM annihilations within subhalos 
can get around these constraints. Specifically, if most of the observed lepton excess comes from 
annihilations in a nearby (within 2 kpc) subhalo along a line of sight toward the galactic center, it 
is possible to match both the lepton and gamma ray observations. We demonstrate that this can 
be achieved in a simple class of particle physics models in which the DM annihilates via a hidden 
leptophilic U(l) vector boson, with explicitly computed Sommerfeld enhancement factors. Gamma 
ray constraints on the main halo annihilations (and CMB constraints from the era of decoupling) 



require the annihilating component of the DM to be subdominant, of order 10 2 
DM density. 
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1. INTRODUCTION 

Nongravitational signals of Dark Matter (DM) have 
been sought after for some time now by the astrophysi- 
cal and particle physics communities. At the same time 
results from the Payload for Antimatter Matter Explo- 
ration and Light-nuclei Astrophysics (PAMELA) experi- 
ment and from the Fermi space telescope suggest a local 
excess positron fraction e + /(e + + e~) at energies above 
10 GeV as well as an excess of e + + e~ peaking around 
500 GeV. Standard cosmic ray propagation models do 
not account for these excesses. An attractive explana- 
tion is that a DM WIMP (weakly interacting massive 
particle) is present in our galaxy at large enough concen- 
trations to self-annihilate into standard model leptons. A 
TeV-scale WIMP annihilating to electron-positron pairs 
could produce such signals. In order to be consistent 
with the observed relic abundance of DM, the annihila- 
tion cross-section (<tv)q ~3x 10 -26 cm 3 s _1 would have 
to be enhanced by a factor of order 100, for example by 
a velocity-dependent Sommerfeld enhancement. 

Many authors [1-10] have explored this possibility, and 
have constrained the allowable mass versus boost factor 
parameter space. However these papers assume that the 
dominant source of indirect signals is from annihilations 
in the main DM halo. In a previous work [11] we consid- 
ered the possibility adding the effects of dark matter sub- 
structure to the theoretical model and we found examples 
where annihilations in subhalos could provide a signifi- 
cant fraction of the observed lepton excesses. We showed 
that one could find a better overall fit to the electron- 
positron data from the Fermi and PAMELA experiments, 
and we suggested that gamma ray constraints which are 
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now putting considerable pressure on these models could 
be alleviated. Our purpose in the present work was to 
ascertain whether this is indeed the case. 

The constraints mentioned come from recent gamma 
ray observations of the galaxy and from Cosmic Mi- 
crowave Background (CMB) measurements. As high 
energy electron-positron pairs are produced and diffuse 
throughout the galaxy, they will emit final-state radia- 
tion as well as scatter on the ambient photon field, giv- 
ing rise to ~ 1-100 GeV gamma rays that should be 
detectable. Given the large expected concentrations of 
both DM and radiation near the galactic center (GC), 
gamma rays from inverse Compton scattering (ICS) near 
the GC are particularly constraining. The Fermi Large 
Area Telescope (LAT) is specifically designed to detect 
gamma rays in this range, and its latest results have been 
used to rule out large regions of parameter space for an- 
nihilating WIMP models [5, 12-14]. 

However in this work we will show that if a sizeable 
proportion of the leptons from DM annihilation originate 
from nearby subhalos, the constraints from GC gamma 
rays can be relieved. Final-state (bremsstrahlung) radi- 
ation from subhalos has been examined by other authors 
[8, 15-20], and ref. [21] has studied the e + + e~ spectrum 
from a nearby subhalo. In this follow-up work we ex- 
tend our previous findings to a prediction of the gamma 
ray spectrum including a full calculation of ICS radia- 
tion in the galaxy, which we compare to the full-sky data 
from the Fermi LAT. We include the expected contribu- 
tion to the gamma ray background coming from back- 
ground electrons and positrons. Using a fully-numerical 
approach, we find that there is less room for new con- 
tributions from the annihilation products of the DM, 
making the constraints on the DM models more severe. 
This is a serious issue even for less cuspy and cored DM 
profiles, that have been shown to satisfy the constraints 
in previous semi-analytic treatments which ignored the 
background gamma ray fluxes. 

In our previous paper we focused on the contributions 
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of distant subhalos to the flux of leptons at Earth. Even 
though these new contributions can improve the fit to the 
lepton data alone, here we show that they do not soften 
the gamma ray constraints sufficiently to be viable. In- 
stead, we focus on the possibility that an accidentally 
nearby subhalo could provide the bulk of the leptonic 
flux. The associated gamma rays would be sufficiently 
hidden by strong backgrounds if this subhalo happened 
to lie between us and the galactic center. The effects 
of nearby subhalos have been previously considered by 
ref. [22] , but only allowing for purely astrophysical boost 
factors, due to the density of the subhalos. Here we find 
that velocity-dependent Sommcrfeld enhancement is cru- 
cial for obtaining a positive outcome. It is precisely be- 
cause of the larger boost factor available within subhalos 
(which have orders of magnitude smaller velocity disper- 
sion) relative to the main halo that we are able to soften 
the gamma ray constraint due to the main halo near the 
GC, yet have a large enough lepton signal from a nearby 
subhalo. In addition, we must assume that the leptophilic 
component of the DM responsible for these processes is 
subdominant to the main inert (for our purposes) com- 
ponent, in order to sufficiently reduce the effective boost 
factor for annihilations in the main halo [14]. This gives 
rise to the interesting possibility that different kinds of 
DM are responsible for the cosmic ray anomalies than 
those which might manifest themselves in direct detec- 
tion experiments. 

Using a modified version of the cosmic ray propaga- 
tion code GALPROP and the data from the recent Via 
Lactea II simulation of dark matter evolution and col- 
lapse in a Milky Way-sized galaxy, we modelled the two- 
dimensional axisymmetric distribution of electrons and 
positrons in the galaxy. These results were combined 
with simulated interstellar radiation field (ISRF) data in 
order to compute a realistic skymap of the gamma ray 
spectrum expected from DM annihilation in the Galaxy, 
which was in turn compared with a year's worth of diffuse 
gamma ray observation from the Fermi LAT. 

We start with a summary the cosmic ray model and re- 
sults of our previous work in Section 2, before discussing 
the relevant ICS and gamma ray physics in Section 3. 
In Section 4 we describe our methodology, and present 
model-independent fits to the data in several scenarios 
for the distribution of subhalos and the halo profiles. In 
particular, we show that an accidentally nearby subhalo 
can provide a promising loophole to the gamma ray con- 
straints on cuspy profiles. We also predict the gamma ray 
flux from the subhalo, which could provide a test of the 
model if future measurements and understanding of back- 
grounds are improved. In section 5 we then demonstrate 
that the boost factors required for this scenario can be 
explicitly realized in a simple class of hidden sector par- 
ticle physics models. We conclude with a discussion of 
the overall viability of this picture in section 6. 



2. COSMIC RAY PROPAGATION 

Inside the galactic diffusion zone, particles and nuclei 
propagate according to the diffusion- loss equation [23], 
which applies to electrons and positrons as follows: 1 

^Ve±(x,p,t) = Q e± (x,£) + V- (D(E) V^(x,p,t)) 

+ ^[&(x,£)V> ± (x,p,t)] . (1) 

?A e ± (x, p, t) denotes the particle number density per unit 
momentum |p|, Q represents the source function, D{E) 
is the spatial diffusion coefficient and 6(x, E) is the en- 
ergy loss coefficient. We seek the steady-state solution of 
equation (1): dip e ± (x, p, t)/dt — 0. 

Since (1) is linear, the leptons from DM annihilation 
travel independently in the astrophysical background. 
The source Q e ± comes from DM annihilation which de- 
pends on the particle physics and the local density of the 
dark matter: 

(2) 

where the prefactor 1/2 is a symmetry factor for self- 
annihilation, udm^E) is the DM energy density, 
(erv)o = 3 x 10~ 26 cm 3 s _1 is the benchmark value for 
standard cosmology to explain the relic density of DM, 
and dA e ± /dE is the energy spectrum of the annihila- 
tion products. Neglecting the effect of soft photons, 
the spectrum can be approximated by the simple form 
dN e ±/dE = 2M D 1 M Q(M DM - E), where Q(x) is the 
usual Heaviside step function, and the factor 2 arises 
because that the final state has two electrons or two 
positrons. The latter has the correct qualitative shape, 
and is easier to implement in GALPROP than would be 
a more exact spectrum. BF denotes the boost factor due 
to Sommcrfeld enhancement, originating from a nonper- 
turbative <~ 1/v correction due to the slow (v/c < a) 
motion of the DM particles. 

To simplify our analysis, we take the boost factor 
BF to be constant throughout the main halo, and tune 
it to provide the best possible fit to available electron 
and positron data. Since the Sommcrfeld effect depends 
strongly on velocity, typical subhalos, which have a much 
smaller velocity dispersion, have a much higher BF, and 
we treat it as an additional free parameter. Although 
each subhalo has different values of BF, we represent 
the subhalo BF by a single average value in this first 



1 The full transport equation also includes the effects of convection 
and diffusive reacceleration, which are mainly important for the 
propagation of heavier species. Here we leave these terms out for 
clarity, although they were included in our full calculations with 
GALPROP. These are important for determining the abundance 
of secondary electrons and positrons, which come from spallation 
and decay of various species. 



3 



part of our analysis, where the BFs are treated as being 
uncorrelated and best fit values are sought. This is not 
a limitation in the case we will eventually focus upon, 
namely domination of the excess lepton signal by a sin- 
gle nearby subhalo. A further complication is that in 
fact BF has a radial dependence within each halo, be- 
cause the velocity dispersion is a function of r, which has 
been fitted by many-body simulations such as Via Lactea 
II [19]. We will take this into account in section 5.1 by 
averaging BF over the phase space of DM in the halos, 
in order to make contact with the results obtained in this 
model-indcpcndcnt part of our analysis. 

The spatial diffusion coefficient can be parametrized as 
follows [24]: 
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Two widely-used approaches exist for solving the diffu- 
sion equation in the Galaxy: semianalytic and fully nu- 
merical. We chose the latter for Galaxy-scale propaga- 
tion, in part because a numerical approach allows for bet- 
ter control over the spatial dependence of the astrophys- 
ical input, such as energy loss due to inverse Compton 
scattering. GALPROP 50. lp [25] is a publicly available 
software package that solves Eq. (1) with an implicit-in- 
time 2D or 3D Crank-Nicholson scheme. In 2D mode, 
it provides a (r, z) map in cylindrical coordinates of the 
number density of each species within the Galactic dif- 
fusion zone. To constrain the diffusion parameters, the 
ratio of measured secondary-to-primary species such as 
B/C or sub-Fe/Fe can be simulated and fit to observa- 
tions. This was done to a very high degree of accu- 
racy in Ref. [24]. We used results from their best fits: 
D = 6.04 x 10 28 cm 2 s _1 (0.19kpc 2 /Myr), and S = 0.41. 

The full energy loss rate is due to synchrotron radiation 
and inverse Compton scattering: 
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a em is the fine structure constant and ub — B 2 /2 is the 
energy density of the galactic magnetic field, for which 
we used the standard parametrization: 



B(r, z) ~ ll^G • exp 



1*1 



lOkpc 2kpc 



(5) 



the energy densities of the three main components 
of the interstellar radiation field (ISRF): CMB radia- 
tion, thermal radiation from dust and starlight, which lie 
mainly in the microwave, infrared and optical regions of 
the electromagnetic spectrum, respectively. GALPROP 
uses position-dependent maps of ISRF compiled by [26] , 
rather than using a constant energy-loss coefficient com- 
puted from a local average. The latter approach (ex- 
plained in section 3 of [27]) is commonly used in the 
semi-analytic model. While it is indeed quite accurate 




Figure 1: Simulated energy density distribution of the inter- 
stellar radiation field (ISRF) within the Milky Way by [26], 
integrated over energies. Top: starlight component. Bottom: 
IR component, from dust. The CMB component is of course 
uniform throughout the galaxy. Color scale is log(density) in 
arbitrary units. 



when dealing with electrons from a smooth Galaxy-wide 
distribution of dark matter, it is an approximation that 
is less precise when considering the propagation into the 
Galaxy of electrons from DM subhalos outside of the dif- 
fusion zone. We will nonetheless make use of the semian- 
alytic method in Section 4.2, when only local propagation 
will be relevant. The position dependence of the ISRF in 
the Galaxy is presented in Figure 1. Further details will 
be discussed in section 3.2. 



2.1. Via Lactea II and GALPROP 

We assumed that the DM was composed of a single 
Dirac fermion x of mass Mjjm annihilating through the 
channel XX —> BB, followed by the decay B — > e + e~, 
where B is some dark sector gauge boson which could 
also be responsible for the Sommerfeld enhancement. We 
considered two astrophysical models for the DM distri- 
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bution: a main halo-only (MH) scenario, in which only 
a large, spherical halo contributed annihilation products; 
and a subhalo (MH+SH) scenario, where the overden- 
sities formed by DM substructure were responsible for 
extra annihilation of DM into electrons and positrons. 
In both cases, we used a spherically symmetric Einasto 
profile for the DM density distribution: 



PEin(0 = p s exp 



1 



(6) 



r is the radial coordinate from the center of the halo, 
p s is the density at r = r s , the distance at which the 
slope p'/p = —2. These parameters are simply related 
to the radius and rotational velocity of a given subhalo 
as explained in Ref. [19]. The shape parameter a can be 
read off from curve-fitting the distributions from TV-body 
simulations such as [28, 29]. It is generally taken to be 
around a ~ 0.17. We took r s — 25 kpc for the main 
galactic halo, with a local dark matter density pq = 0.37 
GeV cm~ 3 in agreement with Via Lactea II and with 
other recent estimates, e.g., [30]. It should be noted that 
many authors use the convention p & = 0.3 GeV cm~ 3 . 
This leads to a factor of (0.3/0.37) 2 = 0.66 diffcrc nee in 
the constraints on the annihilation cross sections, but it 
is of no consequence when it comes to excluding mod- 
els, since constraints come from the ratio of gamma rays 
lepton fluxes, which both scale linearly with pf 3 {av). 

It has been argued that direct observations of rotation 
velocities in the Milky Way are consistent with cored 
DM profiles (see for example ref. [31]). Two such exam- 
ples are the isothermal and Burkcrt [32] ansatzes. The 
Burkert profile has been fitted to the rotation curves of 
galaxies other than our own, but we are not aware of ref- 
erences which attempt to fit the Milky Way. To allow for 
the alternative possibility of a cored main halo, we will 
therefore restrict our attention to the isothermal profile 
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1 + (r/r s y 



(7) 



adopting the values r s = 3.2 kpc and p s = 3.0 GcV/cm 3 
similar to those used by ref. [12]. These values are mo- 
tivated by the constraint on the observed solar density 
Pq (which we take to be somewhat higher than in [12]) 
and on the mass of the Galaxy with 50 kpc as deter- 
mined from circular velocity measurements. However for 
the subhalos we will in all cases assume the Einasto form 
that is suggested by Via Lactea II. 

Via Lactea II [28] was a billion-particle simulation that 
tracked the evolution and collapse of 10 9 particles over 
the history of a Milky Way-sized structure. Data about 
the main galactic halo and the 20,047 largest subhalos 
that the particles (each taken to have mass 4,100 M Q ) 
merged into over the course of the simulation are avail- 
able to the public. While the visible galaxy is only 40 
kpc across, these subhalos extend as far out as 4000 
kpc from the GC. We used the Via Lactea II sub- 
halo data as a model for substructure sourcing electrons 



and positrons (from DM annihilation) at the bound- 
ary of the GALPROP diffusion zone, with an overall 
tunable boost factor for the subhalo annihilation rate. 
In addition to a larger Sommcrfeld enhancement from 
smaller velocity dispersions within each subhalo, we ex- 
pect sub-substructure unresolvable from numerical simu- 
lations to give rise to further enhancement of the annihi- 
lation cross-section. Recent estimates [8] show that such 
sub-subhalos alone could increase annihilation rates by 
as much as a factor of 10. 

Electrons from an extragalactic source have a very par- 
ticular density profile. While the annihilation products 
from the main halo follow a roughly symmetric distribu- 
tion about the GC, SH electrons sourced from the diffu- 
sion zone boundary tend to form a diffuse "shell" near 
the edge of the diffusion zone, as illustrated in Fig. 2. 
Ambient radiation prevents high-energy particles from 
reaching the GC, trapping them near the edge of the 
Galaxy. The large number of subhalos combined with a 
large boost factor can allow some particles to make their 
way to earth, albeit with a fraction of their initial energy. 

We compared the best-fit combination of DM mass and 
boost factor for the MH scenario with the best fits for the 
MH+SH scenario in [11]. The results are summarized in 
table I: a much better fit could be obtained by including 
subhalos and a dark matter particle with Mbm = 2.2 
TeV, rather than the standard MH-only M DM = 1 TcV. 
Of course, the fits are further improved by allowing the 
normalizations of the background electrons and positrons 
to be additional free parameters, denoted as the "freely 
varying background," as opposed to the standard back- 
grounds resulting from GALPROP simulations which in- 
clude the effects of heavier nuclear species. Assuming this 
extra freedom has been advocated or used by numerous 
authors [4, 5, 12, 33]. In table I we also show the fit 
we obtain in the present analysis for the main-halo-only 
case with an isothermal profile and fixed background. It 
is significantly worse than the corresponding one for an 
Einasto profile. 



2.2. Annihilation channels 

While we have mostly focused on the 4c final state, 
there is no reason for other, heavier particles not to be 
produced if the mass of the intermediate gauge boson is 
large enough. Since the amount of Sommerfeld enhance- 
ment ultimately depends on this mass, it is important 
to include the decays to muons and pions. The possible 
final states are all the four-particle combinations of 2e, 
2/i and 2tt. The muon and pion spectra are given by Ref. 
[34] , whose authors were kind enough to provide us with 
the appropriate GALPROP implementation. 

The branching ratios are given by rj = fi/^2fi, where 
the ft are given by 



h - \f, 
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MH+SH 
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14.2 
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Fixed GALPROP background (Einasto) 


MH (4e) 


1.0 


8.2 


144 
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110 




MH+SH 


2.2 


2.1 


175 


177 


146 


1946 


MH (e,/4,7r) 


1.2 


3.8 


109 


112 


118 




Isothermal profile (fixed background) 


MH (4e) 


1.0 


9.1 


186 


195 


113 




MH (e,/4,7r) 


1.2 


3.0 


151 


154 


119 





Table I: First four rows: best fit results from [11], assuming 
Einasto profile. By varying the boost factors of the main halo 
and faraway subhalos separately, we found that the fit to the 
PAMELA and Fermi data from MH annihilations alone could 
be improved by inclusion of SH annihilations as shown. Last 
two rows: new fit for isothermal profile (r s = 3.2 kpc, p s — 
3.0 GeV/cm 3 ), main-halo-only scenario from this work, using 
the fixed GALPROP background, and same parameters as in 
[11]. We assume the annilation to the 4e final state, except 
in the cases "MH (e, n, 7r)" which indicates the the process 
XX BB — > 4£, where £ stands for e , fj, or ir , with 
branching ratios r e = r M = 0.45 and r w — 0.1 as explained in 
Section 2.2. 




5 10 15 20 



r (kpc) 




5 10 15 20 



r (kpc) 



In each the square root factor comes from the phase 
space, while the rest is from the squared matrix ele- 
ment for the decay. Below threshold, fi is defined to be 
zero. For a gauge boson with a mass /i >, 1 GeV, we find 
r e = = 0.45 and r w = 0.1. In this case the electrons 
produced from the final decay of the /x's and 7r's peak at 
a lower energy, thus requiring a slightly higher mass of 
M DM = 1.2 TeV in order to fit the Fermi and PAMELA 
data. This is much smaller than the well-known A/dm — 
2.2 TeV best fit in the pure-muon final state [4, 5, 33] be- 
cause of the large fraction of gauge bosons still decaying 
directly to high-energy electrons. These results are also 
shown in Table I. 



Figure 2: Simulated steady-state distribution of electrons and 
positrons from DM annihilation within the Milky Way diffu- 
sion zone. The galactic center is located at z — 0, r — 0; 
red corresponds to high densities, blue to low densities. Top: 
leptons from the main halo only. Bottom: leptons from the 
subhalos only, sourced from the diffusion zone boundary. Note 
that the scales are different: the peak main halo density (at 
the GC) is about 200 times larger than the peak subhalo den- 
sity (near the edge of the diffusion zone) 



The astrophysical and particle physics dependences of 
each flux can be factorized as 



1 (av) p% dN - 



2 4vr 
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(10) 



3. GAMMA RAY COMPUTATION FROM 
INVERSE COMPTON SCATTERING AND 
BREMSSTRAHLUNG 
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3.1. "Prompt" gamma ray emission 
(bremsstrahlung) 



In each case, the Ji factor depends only upon astrophys- 
ical inputs. The main halo J factor is defined as a line 
of sight (l.o.s.) integral of flux at each pixel: 



Prompt gamma ray emission appears in the final stage 
of DM annihilation, softening the lepton spectrum. The 
flux can be divided into main halo and subhalo parts: 
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In the case of flux originating from many distant subha- 
los, we may treat each one as a point source of radiation. 
In this case, the diffuse flux per solid angle requires a 
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sum over each contributing source with density pi and 
distance di within the observed solid angular region Afl: 



J SI 
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/ 1 



dV 
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(13) 



This clearly depends not only on the density profiles, but 
also on the distribution of subhalos in the Galaxy. We 
will not present the results of the disant subhalo calcula- 
tion of final-state radiation here, since it has been thor- 
oughly explored by other authors in similar contexts. We 
direct the interested reader to references [8, 19, 35]. 

Finally, if a particular subhalo is close enough to sub- 
tend an angle larger than the detector's pixel size, it can 
no longer be treated as a point source: eq. (12) must be 
used, including the angular dependence of the projected 
density profile of the given subhalo, psh (R, </>)■ We 
will return to this case in Section 4.2. 

The particle physics contribution to (10) and (11) 
comes from the photon spectrum, defined as: 
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In the case of a two-lepton final state [36] : 
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where x = 2E~ I / 'y/s and s is the standard Mandelstam 
variable. We are interested in the case of TeV dark mat- 
ter x annihilating to a four-lepton final state, with a 0(1) 
GeV leptophilic gauge boson B as the messenger. The 
annihilation is dominated by \X BB, where the -B's 
are on shell. The cross section can be obtained by first 
computing in the rest frame of the B using the decay 
B — > e + + e~ and then boosting to the lab frame, in 
which the slowly moving DM particles are approximately 
at rest. This can easily be done numerically. We present 
the resulting spectrum in fig. 3. Since we will not make 
use of the final-state bremsstrahlung for other annihila- 
tion channels (4/i or Ait) we will not discuss their spectra. 



3.2. Inverse Compton Scattering 

Charged particles travelling through the interstellar 
medium scatter off ambient photons of the interstellar 
radiation field (ISRF), which is composed of microwave 
(~ 10~ 3 eV) radiation from the cosmic microwave back- 
ground (CMB), infrared (~ 10~ 2 eV) radiation from 
dust, and optical (~ eV) photons from starlight. Along 
with the galactic magnetic fields, this is the main source 
of energy loss for electrons diffusing within the Galaxy. 
We will show that ISRF photons that have scattered 
with TeV-scale electrons have spectra that peak at sev- 
eral hundred GeV, which should fall squarely within the 
measurement window of diffuse gamma rays by the Fermi 
Large Area Telescope (LAT). 




Figure 3: Spectrum of prompt gamma rays (brehmsstrahlung) 
from leptons produced by DM annihilation, as a function of 
x — 2Ej/y/s = E^/Mdm- The red line (upper) represents 
the result of the 2e final state, and the blue line (lower) cor- 
responds to 4e final states. 



Once integrated over scattering angles, the well-known 
Klein-Nishina formula for the Compton scattering pro- 
cess e ± 7 — > e ± 7 / can be integrated along the line of sight 
to give the total flux of scattered photons per solid angle 
arriving on a detector [5, 37]: 
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J, ds represents the line-of-sight integral from the ob- 
server's position to infinity (practically speaking, to the 
edge of the diffusion zone). We have used the definitions: 
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We numerically integrated eq. (16) along the line of 
sight, as well as over the incoming particle energies. 
All the quantities in the integrand are known: we used 
the two-dimensional (r,z) distribution of electrons and 
positrons dn e /dE e from DM annihilations produced with 
GALPROP, as discussed in Section 2. For the ISRF, we 
used a realistic two-dimensional photon energy density 
distribution du-y/dE from [26], which is publicly available 
on the GALPROP website. Both distributions assumed 
cylindrical symmetry around the Galactic axis. For each 
galactic latitude-longitude pair, the line of sight integra- 
tion was performed in a three-dimensional sky from the 
Sun's position to the edge of the diffusion zone which was 
taken to extend to a radius r max = 20 kpc and to a height 
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\z\max = 5 kpc above and below the galactic plane. A 
trapezoidal integration step size of 0.1 kpc was found to 
be numerically converged. The values of dn e /dE e and 
du-t/dE at each step were found in the heliocentric coor- 
dinate system by using a bilinear interpolation scheme. 
On top of the DM annihilation products, we used the 
densities of primary and secondary electrons as well as 
secondary positrons to compute the ICS contribution of 
the background lepton field. This had the effect of fur- 
ther constraining the gamma ray background. 

We performed the integration once per grid point on 
an equally-spaced 20° x 20° latitude-longitude grid of the 
quarter-sky in the ranges 9 — [0,7r/2], (f> — [0, 7r]. This 
was sufficient to reconstruct the entire sky, given the sym- 
metry of the data input. 



3.3. Fermi all-sky diffuse gamma ray measurements 

The Fermi Large Area Telescope (LAT) is a high- 
sensitivity gamma ray instrument capable of detecting 
photons in the ~ 30 MeV to > 300 GeV range. It has 
an effective detector area of <~ 8000 cm 2 , a 2.4 sr field 
of view and can resolve the angle of an incident photon 
to 0.15° at energies above 10 GeV. Data from the first 
year of observation are publicly available from the Fermi 
collaboration. 

We used the all-sky diffuse photon file from the Fermi 
weekly LAT event data webpage [38] . This covered obser- 
vations from mission elapsed time (MET) 239557417 to 
MET 272868753 (seconds), corresponding to 55 weeks of 
observation between August 8 2008 and August 25 2009. 
We processed the photon data with the Fermi LAT sci- 
ence tool software, available from the Fermi Science Sup- 
port Center (FSSC) website. We first removed all events 
with a zenith angle greater than 105° to eliminate Earth 
albedo. The data were further trimmed to keep only the 
photons measured during "good" time intervals. We then 
created an exposure cube from the spacecraft data for the 
corresponding period, to account for effective instrument 
exposure. The data were separated into 0.25° x 0.25° lat- 
itude and longitude bins spanning the entire sky, and into 
16 logarithmically separated bins from 100 MeV to 200 
GeV. Uncertainties were assigned according to Ref. [39] . 
We compared our results to the August-December 2008 
10° < |6| < 20° spectrum presented by the Fermi col- 
laboration [39]. The half-year data agreed exactly, while 
adding the extra 8 months to the full 55-week dataset 
changed the picture only very slightly. We rebinncd the 
data into a 40x40 grid, in correspondance with the ICS 
computation. 

Before proceeding to the results of our numerical anal- 
yses, we should note that many factors contribute to the 
theoretical uncertainty. While we were able to repro- 
duce the results of Simet et al. [24] quite closely, there 
are substantial discrepancies between the results of GAL- 
PROP and other methods of solving the transport equa- 
tion. This lack of agreement is further discussed in [11]. 



There is an additional uncertainty in the injection spec- 
trum of primary electrons, which serve, along with sec- 
ondary electrons and positrons from spallation, as the 
astrophysical background to our results. 



4. EMPIRICAL FITS 

As expected, we found that allowing subhalos to con- 
tribute to the overall flux of DM annihilation products 
reduced the flux of expected gamma rays from the galac- 
tic center, while increasing fluxes at higher galactic lati- 
tudes. The most stringent constraints were from the low- 
longitude regions just above and just below the galactic 
plane, where astrophysical sources of gamma rays are 
less prominent, but the DM distribution is still quite 
dense. Specifically, we used the lower right-hand region 
(-9° < b < -4.5°, 0° < I < 9° in Galactic coordinates) 
which was found to be the most constraining, in agree- 
ment with Ref. [33]. 

After including the ICS from background electrons 
and positrons, we found that the boost factor of a main 
halo 1 TeV DM annihilation process cannot violate the 
bound BF < 25 if the signal is to remain below the 
top Fermi LAT error bars. If we extend the constraint 
to $ 7 < Fermi +2a, this condition is only slightly re- 
laxed to BF < 30. In the case of a 2.2 TeV DM candi- 
date, these bounds become BF < 42 and BF < 52 at 
la and 2<T, respectively. While this agrees qualitatively 
with other works [5, 33], we attribute our more strin- 
gent upper bounds mainly to our higher p Q , as discussed 
in section 2.1, to our inclusion of the ICS contribution 
from background electrons and positrons, but mainly to 
the different method used to solve the diffusion equation 
(!)• 

Using the best fit scenario of Ref. [11], the reduction of 
flux was however not enough to overcome the constraints 
from the Fermi observations. This is illustrated in figure 
4, which shows that the MH+SH scenario still violates 
constraints from the data by as much as 4a. On its own, 
the predicted flux exceeded the data at energies above 
100 GeV by at least 2a, while we expect that additional 
constraints from tt° — > 27 decays should also be large 
in this energy range [40] and push predictions from this 
model even farther outside of the observationally allowed 
region. Allowing the background to freely vary (top sec- 
tion of Table I) made no appreciable difference with re- 
spect to gamma rays, and was not enough to satisfy the 
observational constraints. 

Figure 5 illustrates how the ICS gamma ray flux is 
increased at higher galactic latitudes when subhalos are 
included. It should however be emphasized that the pre- 
dicted fluxes in this region of the sky are still well below 
the level of Fermi observations. 
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4.1. Less cuspy dark matter profiles 

In section 2.1 we mentioned the motivations for con- 
sidering less cuspy DM profiles. Many previous works 
studying the ICS constraints have compared the effects 
of cored versus cuspy DM profiles, noting that the con- 
straints are weaker for cored profiles. To better quantify 
exactly how much cuspiness can be tolerated, it is inter- 
esting to vary the parameters of the Einasto profile that 
control this [14, 41]. In particular, larger values of a and 
r s correspond to less concentrated halos. We ran simu- 
lations of the lepton distribution and gamma ray fluxes 
with slightly different parameters for equation (6) while 
keeping the local density constant at p Q = 0.37 GeV 
cm~ 3 . This is illustrated in fig. 6. Flatter profiles with 
a = 0.20 or 0.25, r s = 30 kpc reduce the gamma ray 
fluxes somewhat, but not enough to bring the predicted 
flux to within the observations in the offending energy 
bins between 10 and 100 GeV. The same is true for the 
isothermal profile, whose corresponding results are shown 
in fig. 7. For both cases, the problem arises because the 
predicted background gamma flux is not far below the 
observed flux in the most constraining bins. This leaves 
very little room for the additional contribution from the 
DM decay products ICS signal. 

Increasing the intermediate gauge boson mass to 1 
GeV, and thus allowing a decay to muons and pions 
according to the branching ratios described in Section 
2.2 does not alleviate the problem. Indeed, the la (2a) 
bounds become BF < 23 (< 28) for an Einasto profile, 
and BF < 63 (< 72) in the isothermal case. These fall 
well short of the required BF = 118 to explain the Fermi 
and PAMELA excesses, as long as the DM mass is in- 
creased to Mom = 1.2 TeV. These results are summa- 
rized in the bottom of Table II. The reason ICS con- 
straints are stronger when muons are included is due to 
the nature of the data. Indeed, the peak of the ICS spec- 
trum lines up with the most constraining data point when 
Mom = 1.2 TeV. This provides a stronger than expected 
constraint, relative to the 4e final state at Mom — 1 TeV. 



4.2. Close subhalo 

The above analyses implicitly assume that no single 
subhalo dominates the lepton signal. But if a subhalo 
happens to be very close (within a kpc) to the solar sys- 
tem, the picture changes significantly, since the electrons 
and positrons from the close subhalo can dominate the 
observed flux, and its gamma ray emissions can come 
from a sizable solid angle in the sky. We treat this case 
separately from the previous subhalo scenario, since a 
larger DM mass is no longer required to produce the ob- 
served lepton signal; rather, the small amount of ICS en- 
ergy loss during propagation from a local subhalo means 
that a 1 TeV-scale DM particle appropriately conforms 
to the Fermi e + + e~ measurements. We concentrate on 
the 4e final state channel, although previous results al- 



Subhalo 


r s (kpc) 


Ps 


log BF 


dmin (pc) 


Vmax (km/s) 


1 


0.01 


69 


4.74 


33.9 


2.9 


2 


0.1 


3.46 


4.34 


95.5 


6.7 


3 


3.2 


0.04 


3.76 


178 


22 


4 


0.9 


1.27 


2.35 


165 


36 


5 


1.1 


2.0 


1.70 


170 


55 


Main halo, 4e channel 


Einasto 


25 


0.048 


1.40 
^ 1.48 




201 - 277 


Isothermal 


3.2 


2.32 


, 1.81 
^ 1.88 




201 - 277 



Main halo, 4e + 4fi + 4ir channel 



Einasto 


25 


0.048 


< 


1.36 
1.45 




201 - 


277 


Isothermal 


3.2 


2.32 


< 


1.80 
1.86 




201 - 


277 



Table II: Upper rows: parameters of each subhalo we exam- 
ined. r s and p s (in units GeV cm -3 ) characterize the halo's 
Einasto profile (with a — 0.17), log BF is the logarithm of 
the necessary boost factor in order to obtain the Fermi lepton 
data entirely from the given subhalo and d m i n is the minimum 
distance (in pc) from our position to such a subhalo along 
the sun-GC axis, with the given boost factor, that would not 
exceed the gamma ray observations. V m ax is the maximum 
circular velocity, which appears in the radial velocity disper- 
sion, fig. 12. Lower rows: similar data for the main halo using 
Einasto or isothermal profiles, but log BF denotes the 1 and 
2a upper limits to satisfy gamma ray constraints. 
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Figure 4: Galactic-center ICS gamma ray flux from the re- 
gion -9° < b < -4.5°, 0° < I < 9° for the MH scenario 
(Mdm = 1 TeV), top black solid line, are reduced in the 
MH+SH scenario (Mdm = 2.2 TeV), middle magenta solid 
line, but not enough to overcome constraints from Fermi LAT 
observations, which are violated by as much as 4a. The pa- 
rameters for the Einasto profile are a — 0.17, r s — 25 kpc. 
The background gamma rays (red solid line) include only ICS 
from background electrons and positrons, but clearly con- 
strain the model even more. Further contributions are ex- 
pected from bremsstrahlung, extragalactic gamma rays and 
7t° decays. The latter may dominate the spectrum at these 
energies and are responsible for the hump shape around 1 
GeV [40]. 
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Figure 5: Mid-latitude ICS gamma rays from the region 42° < 
|6| < 47°, 9° < \i\ < 18°. In this case the MH scenario 
(Mdm = 1 TeV), black solid line, predicts fewer ICS gamma 
rays than the MH+SH scenario (Mdm = 2.2 TeV, magenta 
solid line). At these latitudes constraints are much weaker, 
and neither model is ruled out by the observations. 
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Figure 6: -9° < b < -4.5°, 0° < I < 9° region. Similar to 
previous figures, showing how reducing the cuspiness of the 
Einasto profile (eq. (6)) reduces predicted total gamma ray 
signal (magenta line). Here a = 0.20,0.25 respectively and 
r s — 30 kpc. 



low this to be generalized. The solution depends linearly 
on the spectrum dN e /dE, so that the boost factor re- 
quired to explain the observed lepton excess should scale 
in the same way that it does in the main halo scenario: 
BF( e ^^)/BF ie ~ 118/110, as read from Table I. 

Since GALPROP is not easily adapted in its 2D mode 
to include the effects of a highly localized additional 
source term, we adopt a semi-analytic approach to solve 
the diffusion equation (1) for leptons produced in the 
nearby subhalo. Given that the leptons and gamma rays 



Figure 7: -9° < b < -4.5°, 0° < I < 9° region. Similar 
to previous figures, using the cored isothermal profile with 
r s — 3.2 kpc and p s = 3.0 GeV/cm 3 . 
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Figure 8: Grey regions: scatter plot of p 3 versus r s for subha- 
los in the Via Lactea II simulation. Dots represent the main 
halo (MH) and subhalos given in table II. 



in this scenario would be from a local origin, the spa- 
tial dependence of the interstellar radiation and magnetic 
fields becomes much less important. We used the method 
described in ref. [42] , with the same diffusion parameters 
as presented in section 2 (of the present work), but with 
an energy-loss coefficient parametrized by 



K X > E ) = --IT = — 

at t e 



(19) 



with te — 10 16 s GeV characterizing the local energy loss 
rate. 

We sampled subhalos from the Via Lactea II simula- 
tion to identify examples that could allow for simultane- 
ously fitting the PAMEL A/Fermi lepton fluxes and the 



10 



4 V\. 


gamma rays 


electrons 




\ \V 






\ i \ 


\\ 2 \ 




\ \ 


\\ \ 




\ \ 



1Cf 2 icf 1 10° 10 1 

distance from Earth to SH center (kpc) 



Figure 9: Fluxes of gamma rays and e + e~ from the five 
subhalos presented in Table II. The gamma ray fluxes (curve 
labeled by the number of the corresponding subhalo) are at 
E~f = 137 GeV, whereas the leptons are at an energy of 559 
GeV (the peak of the observed Fermi spectrum). In both 
cases, the amplitude is the predicted flux divided by the ob- 
served flux from the Fermi satellite, such that a value of 10° 
means that the predicted flux is equal to the observed value. 
Boost factors in each case (as given in table II) were fixed to 
allow the Fermi lepton signal to be explained entirely by the 
subhalo. The allowed position of each subhalo with respect 
to earth is therefore the region to the right of each gamma 
ray curve, up to ~ 2 kpc where the lepton flux starts to fall. 

Fermi gamma ray fluxes. Four such examples are labeled 
as SH1-SH4 in table II, and a fifth (SH5) is one that we 
have "engineered" by choosing parameters that are close 
to those of SH4, but with a higher density and hence 
higher circular velocity, dynamically related to each other 
by eq. (13) of [19], 

VL, = fv^G Ps r 2 s (20) 

with fy = 0.897. Due to the higher density, SH5 re- 
quires a lower boost factor to produce the observed lep- 
ton signal, and so it represents a kind of best-case sce- 
nario. The distribution of Via Lactea II subhalos in the 
space of (r s ,p s ) is shown as a scatter plot in fig. 8, and 
the five subhalos of interest are highlighted on this plot. 
They are atypical in the sense of needing a higher-than- 
average central density. A further caveat is that such a 
large r s is unlikely at small distances from the GC due 
to tidal disruption. Indeed, subhalos within the visible 
galaxy in the Via Lactea //simulation were of the order 
r s = 0.05 - 0.85 kpc, falling below the 0.9 - 1.1 kpc 
compatible with the most plausible particle physics sce- 
nario discussed in Section 5. 

Each subhalo was situated along an optimal axis, 
namely that connecting the earth to the GC. Such an ac- 
cidental alignment makes it easier to "hide" the gamma 
rays originating from the subhalo since they are coming 
primarily from the same direction as the GC, where the 



background emissions are strongest. This is also the rea- 
son that the most stringent ICS constraints on the main 
halo arise from the regions 4.5° < \b\ < 9° of galactic 
longitude instead of the most central region. However 
in this case we find that the biggest contribution to the 
emission is from final-state bremsstrahlung rather than 
ICS. The latter is found to produce gamma ray fluxes 
that are 3 orders of magnitude smaller than observed. 
This is consistent with the fact that the main source of 
ICS is IR radiation and starlight, which is concentrated 
far from the vicinity of the solar system. 

Results were then compared to the Fermi lepton and 
gamma ray data in order to establish constraints. The 
strictest gamma constraints were at the largest energy 
data point from the Fermi LAT analysis of E = 162 
GeV, because of the shape of the FSR spectrum, which 
rises steadily until ~ 1 TeV. We used a slightly differ- 
ent region of the sky than in our previous ICS analysis, 
4.5° < \b\ < 9°, 9° < \£\ < 18°, because there were not 
enough good data points in this energy bin at lower lon- 
gitudes to constrain the data. We compared the lepton 
prediction to the Fermi measurements at 559 GeV, where 
the observed e + +e~ spectrum is at a maximum deviation 
from a power law. In both cases we included the addi- 
tional constraints from astrophysical backgrounds com- 
puted by GALPROP and by our ICS routine. 

Results are shown in fig. 9. If the single subhalo is 
allowed to saturate the observed lepton signal, fig. 9 gives 
clear bounds (summarized in Table II) on the proximity 
of each subhalo, providing a minimum distance from the 
solar neighborhhod to such a subhalo. So long as the 
boost factor for the main halo remains sufficiently small, 
this scenario can therefore overcome the ICS constraints 
that restricted the standard MH-only model. 



4.3. Astrophysical prediction and extragalactic 
constraints 

In figure 10 we provide an example of the gamma ray 
flux predicted by the close subhalo scenario, as compared 
to the main halo scenario. The gamma ray flux comes 
predominantly from final state radiation rather than in- 
verse Compton scattering of the annihilation products. 
We chose the energy bin E = 23 GeV, which is the most 
constraining for the main halo case. Although both sce- 
narios converge at high latitudes, low latitude measure- 
ments have already ruled out the main halo scenario, and 
provide a way to constrain the model. With more expo- 
sure and precise removal of point sources, the Fermi LAT 
may provide a diffuse background low enough to rule out 
these predictions. As a further test, census experiments 
such as the upcoming Gaia satellite may provide a precise 
enough map of the local gravitational potential to con- 
firm or rule out the presence of such a DM overdensity 
[43]. Direct measurement of such an overdensity would 
however be difficult: a subhalo such as SH5, located at 
a distance that would not saturate gamma ray bounds, 
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Figure 10: Dependence of predicted gamma ray fluxes on 
galactic latitude b, in the region —9° < £ < 9° at E = 23 
GeV, the most constraining energy bin for the main halo sce- 
nario. Black: main halo scenario (Einasto profile, BF = 110) 
Dashed: subhalo 5, as specified in Table II. Background ICS 
is included in both predictions, but signal is dominated by 
final state radiation. Dots are the Fermi data for that region 
and energy. 



would contribute less than 0.1% of the local DM density. 

From previous works, we infer that extragalactic 
bounds on this scenario are not as strong as the ones 
we have computed above. Bounds from dwarf spheroidal 
galaxies could plausibly be important since the veloc- 
ity dispersions are of the same order as what is required 
for our subhalo enhancement, i.e. ~ 10 — 50 km s — 1 
[44]. However, the most stringent Fermi LAT bounds 
[45] from such galaxies put the upper limit on DM anni- 
hilation into a 2/i final state at around BF = 3000 if only 
final-state radiation is considered, and around 300 if ICS 
bounds are included as well. [46] computed the cosmo- 
logical dark matter annihilation bounds for the same 2/i 
final state scenario, and find that BF larger than 300 is 
excluded at the 90% confidence level. This is using the 
results of the Millennium II structure formation simula- 
tion, and is indeed model-dependent. Extrapolation to 
the 4/x scenario is independent of astrophysics. We can 
therefore take the results of [5, 33] who have construced 
bounds on both channels. They show that FSR bounds 
are consistently an order of magnitude weaker in the 4/i 
case, given the softer photon spectrum in this scenario. 
We can therefore take these extragalactic results to be 
far less constraining than the stringent bounds from the 
center of our own galaxy. 

Finally, we verify that this model does not saturate 
bounds on dipole anisotropy of the cosmic ray e + + e~ 
spectrum. The dipole anisotropy can be defined as 
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Figure 11: Dipole anisotropy 8 of the cosmic ray electron 
and positron flux predicted by SH5 if it saturates the Fermi 
excess. Background cosmic ray electrons and positrons are 
included, and taken to be isotropic. S increases monotonically 
with energy from the red line (60 GeV) to the black line (500 
GeV). 



where C\ is the standard dipole power of the measured 
electron and positron flux in the sky. The Fermi LAT 
collaboration [47] have presented upper bounds on this 
quantity. These range from S < 3 x 10 -3 at E e ~ 60 GeV 
up to S < 9 x 10~ 2 at E e ~ 500 GeV. Given a diffusive 
model, this can be computed [47]: 



S = 



3D(E) \Vn e 



(22) 



(21) 



where D(E) is the diffusion coefficient (3) and n e is the 
density of cosmic ray electrons and positrons, including 
astrophysical backgrounds. Taking the background to be 
isotropic, we computed the dipole anisotropy in the case 
of a single close subhalo producing enough electrons to 
explain the Fermi excess. In every case 5 falls well below 
bounds. Results for SH5 are presented in Figure 11. The 
anisotropy rises monotonically with energy, from 60 GeV 
(red line) to 500 GeV (black line). 



5. PARTICLE PHYSICS REALIZATIONS 

In the previous sections we have identified scenar- 
ios where subhalos could provide the observed excess 
PAMELA and Fermi leptons, from a purely phenomeno- 
logical perspective. In particular, certain values for the 
annihilation cross section boost factors are needed for the 
subhalos, and upper bounds for that of the main halo 
(depending upon assumptions about its density profile) 
were derived. It is interesting to ask whether simple par- 
ticle physics models with boost factors from Sommerfeld 
enhancement can be consistent with these requirements. 

The simplest possibility for model building is dark mat- 
ter that annihilates into light scalar or vector bosons, 
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which subsequently decay into leptons. This class of 
models automatically gives a boost factor to the anni- 
hilation cross section, through multiple exchange of the 
boson, resulting in Sommerfeld enhancement. However 
it is not obvious that one can find models with the de- 
sired boost factors for the subhalos and main halo. One 
constraint that limits our freedom is to not exceed the 
measured density of dark matter. It will turn out that 
our mechanism works most naturally if the DM responsi- 
ble for signals in the galaxy is a subdominant component 
comprising some fraction 1// of the total DM population 
[14], with / > 1. 

We focus on the case of a GeV-scale U(l) vector bo- 
son that kinetically mixes with the photon. Such models 
have the advantage of naturally explaining the coupling 
to light leptons, without producing excess antiprotons 
that would contradict PAMELA observations. Let us 
denote the vector's mass by /i and the coupling by g, 
with a g — g 2 /4tt. If M is the DM mass, then the Som- 
merfeld boost factor is controlled by two dimensionless 
parameters: = /j,/(a g M) and e v = v/(a g c), where v 
is the DM velocity in the center of mass frame. A rea- 
sonably accurate approximation to the exact Sommerfeld 
enhancement is given by the expression [48, 49] 



S = 



sinh X 



€v cosh X — cos 



X 2 



(23) 



where = (7r/12)e^ and X — e v /e^,. (The cosine be- 
comes cosh if the square root becomes imaginary.) 

To take into account leptophilic DM that is only a 
subdominant component of the total DM, suppose that 
a g,th is the value of a g that would give the correct ther- 
mal abundance, which scales like the inverse annihila- 
tion cross section a~ l oc a~ 2 ; then we can parametrize 
a g = -\/7 a g,th- The rate of annihilations goes like 
pfa oc 1/f if pi stands for the leptophilic component of 
the DM. We accordingly define an effective boost factor 



s = 7 



(24) 



where S is the intrinsic Sommerfeld enhancement factor. 
Thus any constraint on S in a theory with / = 1 becomes 
a constraint on S in the more general situation. 



5.1. Averaging of boost factor 

Of course, the DM velocity has no definite value; in- 
stead we need to average over the possible values within 
the subhalos and the main halo, weighted by the appro- 
priate distribution function. We take it to be Maxwell- 
Boltzmann with a cutoff at some escape velocity, 



f(v) = Ae" 3 " 2 / 2 "' 6(v csc - v) 



(25) 



This isotropic form is only an approximation since the 
true distribution has some small anisotropy between the 
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Figure 12: Radial velocity dispersion of subhalos in the Via 
Lactea II simulation, taken from ref. [19] 



radial and angular components; we will for simplicity 
ignore this complication. The velocity dispersion v s — 
(v 2 ) 1 / 2 depends upon the radial distance r from the cen- 
ter of the halo or subhalo. The dependence has been 
measured for the subhalos in the Via Lactea II simula- 
tion; see figure 12. The shape is universal, but is scaled 
along the respective axes by parameters V mayi and ry max 
that depend upon the subhalo. The latter is related to 
the scale radius by ry mo<D — 2.2Vlr s \ the former is given 
by (20) and also listed in table II for the subhalos of 
interest. For numerical purposes we fit the sides of the 
curve passing through the points of fig. 12 by lines (omit- 
ting the rightmost point), and the middle by an inverted 
parabola. 2 We use the same form of v s for the main halo, 
with r s ~ 25 kpc and V" max = 201 km/s. Other authors 
have advocated higher values of the velocity dispersion, 
v s = 309 km/s at r = r Q [50], which would correspond 
to VJnax = 277 km/s in the present parametrization. We 
will also consider the higher value to take account of this 
uncertainty. 

The escape velocity can be computed explicitly 
for the subhalos from the standard result ^v 2 sc = 
G \M(r)/r 2 ) dr, where M(r) = 4tt J^r 2 pdr is the 
mass within radius r. The result for an Einasto profile is 



8tt 



3/c 



7( r (!)- r (M(fi" 



l/a 



( 1 l(JL) a 



(27) 



where r(s,x) is the upper incomplete gamma function. 
For the main halo, this procedure would not be correct 
because of the significant contribution of baryons, not 
included here. We adopt the result for v csc of ref. [14] for 
the main halo (see appendix C of that reference). 



2 The velocity dispersion curve is fit by 

! 1.309 + 0.232a;, 
0.976 - 0.3437x, 
0.9618 - 0.5475a; - 0.4413a; 2 , 

where x = log 10 r/ry mlx and y = v s /V w _ 



x < -0.841, 
x > -0.383 
in between 



(26) 
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Figure 13: Value of gauge coupling leading to correct ther- 
mal relic DM density, a g ^/M , versus squared charge of dark 
Higgs bosons in U(l) model, for several values of DM mass 
M. 



With these ingredients, we can compute an average 
Sommerfeld enhancement factor (S) for each subhalo: 



(S) 



C drr 2 p 2 / <f Vl d?v 2 f{v x ) /(« 2 )^(||«i - v 2 \) 



f 

J 1 



r2 dr r 2 p 2 



(28) 

The factor of | in the argument of S occurs because the 
v appearing in eq. (23) through e v is half of the relative 
velocity, p 2 is the appropriate weighting factor because 
the rate of annihilations is proportional to (av)p 2 . For 
the subhalos, the range of integration for r is from 
to oo, but for the main halo we take lower and upper 
limits ri.2 that correspond to the angular region of the 
sky that is used to set the gamma ray constraints: r\ = 
0.67 kpc and r 2 = 1.34 kpc. The reason is that the 
bound S < 30 for the main halo comes from the gamma 
ray constraint rather than from lepton production. We 
are thus interested in the boost factor relevant to the 
region 4.5° < \b\ < 9° of galactic latitude. The distances 
of closest approach to the galactic center, hence largest 
rate of 7 ray production associated with these lines of 
sight, are given by r = r sin b. 



5.2. Relic Density Constraint 

The enhancement factor (23) depends rather strongly 
on the gauge coupling a g ; therefore it is interesting to 
know what constraint the relic density places upon a g . 
The effect of a Sommerfeld-enhanced DM model on the 
relic densitie has been discussed by [51]. Notice that DM 
transforming under a U(l) gauge symmetry as we have 
assumed must be Dirac and therefore could have a relic 
density through its asymmetry, similar to baryons. How- 
ever, unless the DM was never in thermal equilibrium, 
then a g should not be less than the usual value a Qt th 
leading to the correct relic density, since otherwise the 
thermal component will be too large. 



There are two kinds of final states for annihilation of 
DM in this class of models: into a pair of gauge bosons 
B^, by virtual DM exchange in the t and u channels, or 
into dark Higgs bosons h, by exchange of a gauge boson 
in the s channel. Assuming the DM (x) is much heavier 
than the final states, the respective squared amplitudes, 
averaged over initial and summed over final spins, are 



1 -T\M\ 2 =h^ i+2v2 ^ 



XX -> BB 
XX hh 



(29) 

where q is the U(l) charge of h relative to x (replace 
q 2 —> Ylili f° r multiple Higgs bosons), 8 is the scat- 
tering angle, and we have included the leading depen- 
dence on the initial velocity v in the center of mass 
frame. The factor cos 2 8 averages to 2/3 in the inte- 
gral over 8. In computing the associated cross section, 
it must be remembered that the 2B final state consists 
of identical particles, while the Higgs channel does not. 
The total amplitude can therefore be written in the form 
i£|.M| 2 = 4g 4 {a + bv 2 ), with 



a = l + |Ei«?, & = 2(1- &£i«? 



(30) 



if we use the phase space for identical particles. 

To find the cross section relevant during freeze-out in 
the early universe, we thermally average the w-dependent 
<7i> re i following rcf. [52]. We include approximately the 
effect of Sommerfeld enhancement as described there, to 
obtain 



2M 2 
T 
M 



(31) 



The terms that are subleading in a gi but enhanced by 
a/ M/T, are due to the Sommerfeld correction. We ap- 
proximate the freezeout temperature as T = M/20, the 
usual result of solving the Boltzmann equation for DM 
in the TeV mass range, and equate (crv le i) to the value 
(<jv)o = 3 x 10~ 26 cm 3 /s usually assumed to give the 
correct relic density. This gives an implicit equation for 
oi g ,tb of the form a 2 = ciM 2 (av)o/(l + c 2 a g ), which how- 
ever quickly converges by numerically iterating. Fig. 13 
displays the resulting dependence of a g ,th/M on Y^iQi 
for several values of M . 

The bound that the density of the leptophilic DM com- 

g,th- 



ponent not exceed the total DM density is a„ > a 



We parametrize the coupling by a g = y/J a 9i th with 
/ > 1 in what follows. 



5.3. Interpolation between 4e and mixed final states 

In our numerical computations with GALPROP, we 
considered two cases for the final state annihilation chan- 
nels: either XX ~^ 4e, applicable for gauge bosons with 
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Figure 14: Solid lines: predicted main halo boost factor for 
thermal value of a g , with dark Higgs boson charges ]TV qf = 
16 and maximum circular velocity Knax = 277 km/s. Up- 
per curve is for Einasto profile, lower for isothermal. Dashed 
line is 2<r upper limit from gamma rays produced by inverse 
Compton scattering. The failure to satisfy this bound even 
with large dark Higgs content and large Vmax drives us to 
consider larger than thermal gauge couplings, / > 1. 

mass fj, < 2m M , or to a mixture of electrons, muons and 
charged pions, appropriate for decays of gauge bosons 
with mass greater than 2m 7r . The relative abundance 
of e, fi and ir in the mixed final state can be computed 
from the branching fractions of the decays, discussed in 
connection with eq. (8). 

For intermediate values of the gauge boson mass, 
2m M < fi < 2m 3r , we can use the branching ratios to in- 
terpolate between our maximum- allowed MH or best-fit 
SH boost factors for the 4e case and those of the fidu- 
cial e + /i + 7T case. The maximum allowed boost factors 
of the main halo complying with the ICS constraints are 
taken from table II. To estimate the best fit boost factors 
for the subhalos in the fiducial e + fi + n final state, we 
rescale the 4e results shown in table II by the ratio of 
best-fit boost factors for the main halo, in the MH-only 
scenario. These ratios are 118/110 for the Einasto pro- 
file and 119/113 for the isothermal, quite close to unity, 
and so the best-fit values of the SH boost factors hardly 
depend upon this scaling. More significant is the change 
in the best-fit mass, from M = 1.0 to 1.2 TeV, which 
enters into the computation of the Sommerfeld enhance- 
ment and the value of the gauge coupling (a g ~ M). We 
use the branching ratios to interpolate M as well. For 
the MH upper bounds in the small- and large-/z regions, 
we use the values from Table II, and interpolate similarly 
for intermediate \i. 



5.4. Theoretical fits 

For a given value of the gauge coupling a g , we can de- 
termine the predicted boost factors as well as the desired 
values for each subhalo, as a function of the gauge boson 



mass fi, and similarly for the main halo, except here we 
have an upper bound on (S) rather than a best-fit value. 
This bound in fact presents the biggest challenge to find- 
ing a working particle physics model. For a g close to 
the thermal relic density value a g ^h, the predicted boost 
factor of the main halo far exceeds the bound (S) 30, 
even if we try to decrease (S) by reducing a g via a large 
hidden Higgs content or by increasing the dispersion of 
the main halo. Fig. 14 illustrates the discrepancy for 
J2 i q? = 16 and V ma x = 277 km/s. Lower values of V ma x 
or '^2 i qf only make this tension worse. 

As we mentioned above, even though it is not theoreti- 
cally possible to make the gauge coupling weak enough to 
solve this problem, ironically one can rescue the scenario 
by increasing a g beyond the thermal value, since this 
suppresses the relic density of the DM component we are 
interested in, and thus reduces the scattering rate. Al- 
lowing a g — V/c^g.th decreases both the density of the 
leptophilic component and the effective boost factor by 
!//• With / ~ 50 — 500, depending upon the shape of 
the main halo DM density profile, we can satisfy the con- 
straint on the MH and still have a large enough boost in 
certain hypothetical nearby subhalos for them to supply 
the observed lepton excess. The minimum value of / that 
is needed is larger for a cuspy main halo. 

We give two working examples in figure 15, one with 
/ = 500 and V max — 277 km/s (the larger value advo- 
cated in ref. [50] ) and assuming an Einasto profile for the 
main halo, and the other having / = 50 and Vnax = 201 
km/s (the more standard assumption for the velocity dis- 
persion), with an isothermal halo. In these figures the 
averaged boost factor (S) of the relevant subhalos are 
plotted as solid lines, while the required values of (S) 
are the dashed curves. Wherever these intersect repre- 
sents a possible value of the gauge boson mass to consis- 
tently explain the observed lepton excess. At the same 
time, the main halo boost factor (lowest solid curve in 
the small- fi region) must lie below the black dashed lines 
to satifsy gamma ray constraints. The rationale for tak- 
ing the larger value of Vmax for the Einasto profile is that 
larger velocities help to suppress the boost factors and 
thus make it easier to satisfy the ICS constraint, so that 
we are not forced to choose an even larger value of /. 
The isothermal profile is less constrained. 

In the first panel of fig. 15 with the Einasto profile, 
only subhalos SH4 and SH5 have large enough boost fac- 
tors to ever reach the required values. There are many 
points of intersection, but mainly those for SH5 and in 
the mass range n < 750 MeV are consistent with the 
gamma ray bounds on the main halo. For the isothermal 
halo, these constraints are less stringent, and it is possible 
to find points of intersection using / = 50 for all five of 
the sample subhalos, although they are much more rare 
for SH1-SH3 than for SH4 and SH5. In this example, 
the intersection points that respect the ICS bound are 
restricted to fi < 1 GeV. For larger values of /, all the 
boost factors will be further suppressed, and fi > 1 GeV 
will become allowed for SH4 and SH5. 
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Figure 15: Predicted effective boost factors {§) as a function of gauge boson mass (solid curves) and target values (or upper limit 
in case of main halo, dashed curves) to explain PAMELA/Fermi lepton observations and Fermi gamma ray constraints. Pair 
of dashed curves for main halo (MH) correspond respectively to 1 and 2a upper limits. Left panel is for / = 100, V ma x = 277 
km/s, Einasto main halo profile; right is for / = 25, V^max = 201 km/s, isothermal main halo profile. Subhalos are those of 
table II. Points which satisfy all constraints are those where subhalo curves intersect their corresponding dashed line while the 
main halo curve lies below its dashed lines. 



One advantage of requiring large / is that the corre- 
sponding dilution of the DM density by 1// insures that 
the model satisfies stringent CMB constraints from anni- 
hilations in the early universe changing the optical depth 
[53-55], as pointed out in [14]. The CMB constraint is 
shown in fig. 16, along with the PAMELA/Fermi allowed 
regions from ref. [33] for 4e and final states. The 4e 
case is allowed by the CMB constraint, but 4/i is ruled 
out. Because our model has at most a fraction of 0.45 of 
muons in the final state, it is probably already safe, but 
the additional weakening of the bound by the factor 1// 
ensures that this will be the case. Similarly, our scenario 
overcomes the no-go result of ref. [10], which pointed out 
that Sommerfeld enhanced annihilation in the early uni- 
verse leads to constraints on the MH boost factor which 
are lower than those needed to explain the lepton anoma- 
lies. Our MH boost factor can satisfy these constraints 
since the MH is no longer considered to be the source of 
the excess leptons. 




_25 1 1 1 1 1 1 1 1 1 

2 2.5 3 3.5 4 
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Figure 16: Allowed regions for PAMELA and Fermi excess 
leptons, and upper bounds from inverse Compton gamma 
rays, from ref. [33], for Einasto profile with a = 0.17 and 
r s — 20 kpc. CMB constraint is from ref. [53]. 



6. DISCUSSION AND CONCLUSIONS 



The Sommerfeld enhancement is nearly saturated for 
the low velocities of the subhalos at these large values of 
otg r~j 0.1 — 0.35, so their (S) curves are nearly overlapping 
except at the smallest gauge boson masses. The main 
halo boost factor is not saturated on the other hand, and 
lies below the FSR bound for most values of fx. We have 
chosen the gauge couplings, parametrized by /, to nearly 
saturate the FSR bound. By taking larger a g (larger /), 
the bounds could be satisfied by a larger margin. But 
this would also reduce the (S) values of the subhalos 
by a similar amount, making it more difficult to get a 
large enough lepton signal from SHI— SH3. SH4 and SH5 
would remain robust possible explanations. 



We have shown that gamma ray constraints on 
leptophilic annihilating dark matter arc significantly 
stronger than in previous studies, when we take into 
account the contributions to inverse Compton scatter- 
ing from primary and secondary electrons and positrons, 
before including excess leptons from the DM annihila- 
tion. We attribute part of this difference to the method 
of solving the diffusion equation (1) — fully numeri- 
cal rather than semi- analytic — meaning that the (r, z) 
space-dependence of the diffusion coefficient is taken into 
account. The difference between the predicted and ob- 
served spectra of gamma rays is greatly reduced, leaving 
little room for new contributions. Because of this, even 
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cored halos, which were allowed by other analyses, be- 
come excluded. However, we find that these constraints 
can be weakened and possibly overcome if annihilations 
in a nearby subhalo are the dominant source of anoma- 
lous leptons, rather than annihilations in the main galac- 
tic DM halo. In this way, the PAMEL A/Fermi cosmic 
ray excesses can be explained, without violating bounds 
from the recent Fermi LAT diffuse gamma ray survey. 

It must be admitted that the subhalo loophole we 
present is rather special. First, only atypically dense sub- 
halos, relative to the sample provided by Via Lactea II, 
give a large enough boost factor (see fig. 8). Second, the 
subhalo would need to accidentally line up nearly with 
the galactic center in order for the ICS gamma rays asso- 
ciated with these leptons to be sufficiently hidden by the 
noisy background of the GC. Of course, had we neglected 
ICS contributions of background electrons, similar to pre- 
vious studies, less fine tuning of the subhalo properties 
would be necessary. Also we do not require the subhalo 
to be particularly close; fig. 9 shows that the lepton flux 
only starts to fall at distances of ~ 3 kpc. Our find- 
ing could be regarded as a proof of concept. It is possi- 
ble that the effects of unresolved substructure within the 
subhalo [15], which can increase the boost factor, would 
also make the scenario work more easily. On the positive 
side, there is the opportunity of testing whether there is 
such a nearby subhalo, since we predict the spectrum of 
ICS gamma rays it contributes (see fig. 10). A better 
understanding of backgrounds, for example from point 
sources, could make it possible to rule out the proposal. 

On the particle physics side, we have shown in detail 
that the subhalo scenario can be made consistent with 
one of the simplest models of leptophilic dark matter, 



where the DM is in a hidden sector that communicates 
with the standard model only through kinetic mixing 
with hypercharge of a new gauge boson in the GeV mass 
range. The relative couplings to leptons and charged pi- 
ons are completely specified and the model has only two 
free parameters, the gauge coupling a g and gauge bo- 
son mass n (the DM mass M is fixed by fitting to the 
spectrum of anomalous e + + e~). The gauge coupling is 
constrained by the relic density of the DM. The Sommcr- 
feld enhancement factor is completely fixed by (a g ,fi, M) 
and the kinematical halo properties. We find (similarly 
to ref. [14]) that the predicted boost factor for the main 
halo is always too large to satisfy ICS constraints unless 
the leptophilic component of the DM is small, compris- 
ing a fraction of order 1/f = 0.02 - 0.002 of the total 
DM. The small fraction can be achieved by assuming a g 
is larger than the value required for the usual thermal 
abundance by the factor / <~ 50 — 500. This raises the 
interesting possibility that the DM that may be respon- 
sible for the cosmic ray anomalies is distinct from the 
dominant DM species that might be discovered by direct 
detection. 
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